home *** CD-ROM | disk | FTP | other *** search
/ Aminet 40 / Aminet 40 (2000)(Schatztruhe)[!][Dec 2000].iso / Aminet / dev / lang / Python16_Src.lha / Python16_Source / Objects / complexobject.c < prev    next >
Encoding:
C/C++ Source or Header  |  2000-09-10  |  12.3 KB  |  624 lines
  1. /* Complex object implementation */
  2.  
  3. /* Borrows heavily from floatobject.c */
  4.  
  5. /* Submitted by Jim Hugunin */
  6.  
  7. #ifndef WITHOUT_COMPLEX
  8.  
  9. #include "Python.h"
  10. #include "mymath.h"
  11. #include "protos/complexobject.h"
  12.  
  13. #ifdef HAVE_LIMITS_H
  14. #include <limits.h>
  15. #endif
  16.  
  17.  
  18. /* elementary operations on complex numbers */
  19.  
  20. static Py_complex c_1 = {1., 0.};
  21.  
  22. Py_complex c_sum(a,b)
  23.     Py_complex a,b;
  24. {
  25.     Py_complex r;
  26.     r.real = a.real + b.real;
  27.     r.imag = a.imag + b.imag;
  28.     return r;
  29. }
  30.  
  31. Py_complex c_diff(a,b)
  32.     Py_complex a,b;
  33. {
  34.     Py_complex r;
  35.     r.real = a.real - b.real;
  36.     r.imag = a.imag - b.imag;
  37.     return r;
  38. }
  39.  
  40. Py_complex c_neg(a)
  41.     Py_complex a;
  42. {
  43.     Py_complex r;
  44.     r.real = -a.real;
  45.     r.imag = -a.imag;
  46.     return r;
  47. }
  48.  
  49. Py_complex c_prod(a,b)
  50.     Py_complex a,b;
  51. {
  52.     Py_complex r;
  53.     r.real = a.real*b.real - a.imag*b.imag;
  54.     r.imag = a.real*b.imag + a.imag*b.real;
  55.     return r;
  56. }
  57.  
  58. Py_complex c_quot(a,b)
  59.     Py_complex a,b;
  60. {
  61.     Py_complex r;
  62.     double d = b.real*b.real + b.imag*b.imag;
  63.     if (d == 0.)
  64.         errno = EDOM;
  65.     r.real = (a.real*b.real + a.imag*b.imag)/d;
  66.     r.imag = (a.imag*b.real - a.real*b.imag)/d;
  67.     return r;
  68. }
  69.  
  70. Py_complex c_pow(a,b)
  71.     Py_complex a,b;
  72. {
  73.     Py_complex r;
  74.     double vabs,len,at,phase;
  75.     if (b.real == 0. && b.imag == 0.) {
  76.         r.real = 1.;
  77.         r.imag = 0.;
  78.     }
  79.     else if (a.real == 0. && a.imag == 0.) {
  80.         if (b.imag != 0. || b.real < 0.)
  81.             errno = ERANGE;
  82.         r.real = 0.;
  83.         r.imag = 0.;
  84.     }
  85.     else {
  86.         vabs = hypot(a.real,a.imag);
  87.         len = pow(vabs,b.real);
  88.         at = atan2(a.imag, a.real);
  89.         phase = at*b.real;
  90.         if (b.imag != 0.0) {
  91.             len /= exp(at*b.imag);
  92.             phase += b.imag*log(vabs);
  93.         }
  94.         r.real = len*cos(phase);
  95.         r.imag = len*sin(phase);
  96.     }
  97.     return r;
  98. }
  99.  
  100. static Py_complex c_powu(x, n)
  101.     Py_complex x;
  102.     long n;
  103. {
  104.     Py_complex r, p;
  105.     long mask = 1;
  106.     r = c_1;
  107.     p = x;
  108.     while (mask > 0 && n >= mask) {
  109.         if (n & mask)
  110.             r = c_prod(r,p);
  111.         mask <<= 1;
  112.         p = c_prod(p,p);
  113.     }
  114.     return r;
  115. }
  116.  
  117. static Py_complex c_powi(x, n)
  118.     Py_complex x;
  119.     long n;
  120. {
  121.     Py_complex cn;
  122.  
  123.     if (n > 100 || n < -100) {
  124.         cn.real = (double) n;
  125.         cn.imag = 0.;
  126.         return c_pow(x,cn);
  127.     }
  128.     else if (n > 0)
  129.         return c_powu(x,n);
  130.     else
  131.         return c_quot(c_1,c_powu(x,-n));
  132.  
  133. }
  134.  
  135. PyObject *
  136. PyComplex_FromCComplex(cval)
  137.     Py_complex cval;
  138. {
  139.     register PyComplexObject *op;
  140.  
  141.     /* PyObject_New is inlined */
  142.     op = (PyComplexObject *) PyObject_MALLOC(sizeof(PyComplexObject));
  143.     if (op == NULL)
  144.         return PyErr_NoMemory();
  145.     PyObject_INIT(op, &PyComplex_Type);
  146.     op->cval = cval;
  147.     return (PyObject *) op;
  148. }
  149.  
  150. PyObject *
  151. PyComplex_FromDoubles(real, imag)
  152.     double real, imag;
  153. {
  154.     Py_complex c;
  155.     c.real = real;
  156.     c.imag = imag;
  157.     return PyComplex_FromCComplex(c);
  158. }
  159.  
  160. double
  161. PyComplex_RealAsDouble(op)
  162.     PyObject *op;
  163. {
  164.     if (PyComplex_Check(op)) {
  165.         return ((PyComplexObject *)op)->cval.real;
  166.     } else {
  167.         return PyFloat_AsDouble(op);
  168.     }
  169. }
  170.  
  171. double
  172. PyComplex_ImagAsDouble(op)
  173.     PyObject *op;
  174. {
  175.     if (PyComplex_Check(op)) {
  176.         return ((PyComplexObject *)op)->cval.imag;
  177.     } else {
  178.         return 0.0;
  179.     }
  180. }
  181.  
  182. Py_complex
  183. PyComplex_AsCComplex(op)
  184.     PyObject *op;
  185. {
  186.     Py_complex cv;
  187.     if (PyComplex_Check(op)) {
  188.         return ((PyComplexObject *)op)->cval;
  189.     } else {
  190.         cv.real = PyFloat_AsDouble(op);
  191.         cv.imag = 0.;
  192.         return cv;
  193.     }   
  194. }
  195.  
  196. static void
  197. complex_dealloc(op)
  198.     PyObject *op;
  199. {
  200.     PyObject_DEL(op);
  201. }
  202.  
  203.  
  204. static void
  205. complex_buf_repr(buf, v)
  206.     char *buf;
  207.     PyComplexObject *v;
  208. {
  209.     if (v->cval.real == 0.)
  210.         sprintf(buf, "%.12gj", v->cval.imag);
  211.     else
  212.         sprintf(buf, "(%.12g%+.12gj)", v->cval.real, v->cval.imag);
  213. }
  214.  
  215. static int
  216. complex_print(v, fp, flags)
  217.     PyComplexObject *v;
  218.     FILE *fp;
  219.     int flags; /* Not used but required by interface */
  220. {
  221.     char buf[100];
  222.     complex_buf_repr(buf, v);
  223.     fputs(buf, fp);
  224.     return 0;
  225. }
  226.  
  227. static PyObject *
  228. complex_repr(v)
  229.     PyComplexObject *v;
  230. {
  231.     char buf[100];
  232.     complex_buf_repr(buf, v);
  233.     return PyString_FromString(buf);
  234. }
  235.  
  236. static int
  237. complex_compare(v, w)
  238.     PyComplexObject *v, *w;
  239. {
  240. /* Note: "greater" and "smaller" have no meaning for complex numbers,
  241.    but Python requires that they be defined nevertheless. */
  242.     Py_complex i, j;
  243.     i = v->cval;
  244.     j = w->cval;
  245.     if (i.real == j.real && i.imag == j.imag)
  246.        return 0;
  247.     else if (i.real != j.real)
  248.        return (i.real < j.real) ? -1 : 1;
  249.     else
  250.        return (i.imag < j.imag) ? -1 : 1;
  251. }
  252.  
  253. static long
  254. complex_hash(v)
  255.     PyComplexObject *v;
  256. {
  257.     double intpart, fractpart;
  258.     int expo;
  259.     long hipart, x;
  260.     /* This is designed so that Python numbers with the same
  261.        value hash to the same value, otherwise comparisons
  262.        of mapping keys will turn out weird */
  263.  
  264. #ifdef MPW /* MPW C modf expects pointer to extended as second argument */
  265. {
  266.     extended e;
  267.     fractpart = modf(v->cval.real, &e);
  268.     intpart = e;
  269. }
  270. #else
  271.     fractpart = modf(v->cval.real, &intpart);
  272. #endif
  273.  
  274.     if (fractpart == 0.0 && v->cval.imag == 0.0) {
  275.         if (intpart > 0x7fffffffL || -intpart > 0x7fffffffL) {
  276.             /* Convert to long int and use its hash... */
  277.             PyObject *w = PyLong_FromDouble(v->cval.real);
  278.             if (w == NULL)
  279.                 return -1;
  280.             x = PyObject_Hash(w);
  281.             Py_DECREF(w);
  282.             return x;
  283.         }
  284.         x = (long)intpart;
  285.     }
  286.     else {
  287.         fractpart = frexp(fractpart, &expo);
  288.         fractpart = fractpart * 2147483648.0; /* 2**31 */
  289.         hipart = (long)fractpart; /* Take the top 32 bits */
  290.         fractpart = (fractpart - (double)hipart) * 2147483648.0;
  291.                         /* Get the next 32 bits */
  292.         x = hipart + (long)fractpart + (long)intpart + (expo << 15);
  293.                         /* Combine everything */
  294.  
  295.         if (v->cval.imag != 0.0) { /* Hash the imaginary part */
  296.             /* XXX Note that this hashes complex(x, y)
  297.                to the same value as complex(y, x).
  298.                Still better than it used to be :-) */
  299. #ifdef MPW
  300.             {
  301.                 extended e;
  302.                 fractpart = modf(v->cval.imag, &e);
  303.                 intpart = e;
  304.             }
  305. #else
  306.             fractpart = modf(v->cval.imag, &intpart);
  307. #endif
  308.             fractpart = frexp(fractpart, &expo);
  309.             fractpart = fractpart * 2147483648.0; /* 2**31 */
  310.             hipart = (long)fractpart; /* Take the top 32 bits */
  311.             fractpart =
  312.                 (fractpart - (double)hipart) * 2147483648.0;
  313.                         /* Get the next 32 bits */
  314.             x ^= hipart + (long)fractpart +
  315.                 (long)intpart + (expo << 15);
  316.                         /* Combine everything */
  317.         }
  318.     }
  319.     if (x == -1)
  320.         x = -2;
  321.     return x;
  322. }
  323.  
  324. static PyObject *
  325. complex_add(v, w)
  326.     PyComplexObject *v;
  327.     PyComplexObject *w;
  328. {
  329.     Py_complex result;
  330.     PyFPE_START_PROTECT("complex_add", return 0)
  331.     result = c_sum(v->cval,w->cval);
  332.     PyFPE_END_PROTECT(result)
  333.     return PyComplex_FromCComplex(result);
  334. }
  335.  
  336. static PyObject *
  337. complex_sub(v, w)
  338.     PyComplexObject *v;
  339.     PyComplexObject *w;
  340. {
  341.     Py_complex result;
  342.     PyFPE_START_PROTECT("complex_sub", return 0)
  343.     result = c_diff(v->cval,w->cval);
  344.     PyFPE_END_PROTECT(result)
  345.     return PyComplex_FromCComplex(result);
  346. }
  347.  
  348. static PyObject *
  349. complex_mul(v, w)
  350.     PyComplexObject *v;
  351.     PyComplexObject *w;
  352. {
  353.     Py_complex result;
  354.     PyFPE_START_PROTECT("complex_mul", return 0)
  355.     result = c_prod(v->cval,w->cval);
  356.     PyFPE_END_PROTECT(result)
  357.     return PyComplex_FromCComplex(result);
  358. }
  359.  
  360. static PyObject *
  361. complex_div(v, w)
  362.     PyComplexObject *v;
  363.     PyComplexObject *w;
  364. {
  365.     Py_complex quot;
  366.     PyFPE_START_PROTECT("complex_div", return 0)
  367.     errno = 0;
  368.     quot = c_quot(v->cval,w->cval);
  369.     PyFPE_END_PROTECT(quot)
  370.     if (errno == EDOM) {
  371.         PyErr_SetString(PyExc_ZeroDivisionError, "complex division");
  372.         return NULL;
  373.     }
  374.     return PyComplex_FromCComplex(quot);
  375. }
  376.  
  377. static PyObject *
  378. complex_remainder(v, w)
  379.     PyComplexObject *v;
  380.     PyComplexObject *w;
  381. {
  382.         Py_complex div, mod;
  383.     errno = 0;
  384.     div = c_quot(v->cval,w->cval); /* The raw divisor value. */
  385.     if (errno == EDOM) {
  386.         PyErr_SetString(PyExc_ZeroDivisionError, "complex remainder");
  387.         return NULL;
  388.     }
  389.     div.real = floor(div.real); /* Use the floor of the real part. */
  390.     div.imag = 0.0;
  391.     mod = c_diff(v->cval, c_prod(w->cval, div));
  392.  
  393.     return PyComplex_FromCComplex(mod);
  394. }
  395.  
  396.  
  397. static PyObject *
  398. complex_divmod(v, w)
  399.     PyComplexObject *v;
  400.     PyComplexObject *w;
  401. {
  402.         Py_complex div, mod;
  403.     PyObject *d, *m, *z;
  404.     errno = 0;
  405.     div = c_quot(v->cval,w->cval); /* The raw divisor value. */
  406.     if (errno == EDOM) {
  407.         PyErr_SetString(PyExc_ZeroDivisionError, "complex divmod()");
  408.         return NULL;
  409.     }
  410.     div.real = floor(div.real); /* Use the floor of the real part. */
  411.     div.imag = 0.0;
  412.     mod = c_diff(v->cval, c_prod(w->cval, div));
  413.     d = PyComplex_FromCComplex(div);
  414.     m = PyComplex_FromCComplex(mod);
  415.     z = Py_BuildValue("(OO)", d, m);
  416.     Py_XDECREF(d);
  417.     Py_XDECREF(m);
  418.     return z;
  419. }
  420.  
  421. static PyObject *
  422. complex_pow(v, w, z)
  423.     PyComplexObject *v;
  424.     PyObject *w;
  425.     PyComplexObject *z;
  426. {
  427.     Py_complex p;
  428.     Py_complex exponent;
  429.     long int_exponent;
  430.  
  431.      if ((PyObject *)z!=Py_None) {
  432.         PyErr_SetString(PyExc_ValueError, "complex modulo");
  433.         return NULL;
  434.     }
  435.  
  436.     PyFPE_START_PROTECT("complex_pow", return 0)
  437.     errno = 0;
  438.     exponent = ((PyComplexObject*)w)->cval;
  439.     int_exponent = (long)exponent.real;
  440.     if (exponent.imag == 0. && exponent.real == int_exponent)
  441.         p = c_powi(v->cval,int_exponent);
  442.     else
  443.         p = c_pow(v->cval,exponent);
  444.  
  445.     PyFPE_END_PROTECT(p)
  446.     if (errno == ERANGE) {
  447.         PyErr_SetString(PyExc_ValueError,
  448.                 "0.0 to a negative or complex power");
  449.         return NULL;
  450.     }
  451.  
  452.     return PyComplex_FromCComplex(p);
  453. }
  454.  
  455. static PyObject *
  456. complex_neg(v)
  457.     PyComplexObject *v;
  458. {
  459.     Py_complex neg;
  460.     neg.real = -v->cval.real;
  461.     neg.imag = -v->cval.imag;
  462.     return PyComplex_FromCComplex(neg);
  463. }
  464.  
  465. static PyObject *
  466. complex_pos(v)
  467.     PyComplexObject *v;
  468. {
  469.     Py_INCREF(v);
  470.     return (PyObject *)v;
  471. }
  472.  
  473. static PyObject *
  474. complex_abs(v)
  475.     PyComplexObject *v;
  476. {
  477.     double result;
  478.     PyFPE_START_PROTECT("complex_abs", return 0)
  479.     result = hypot(v->cval.real,v->cval.imag);
  480.     PyFPE_END_PROTECT(result)
  481.     return PyFloat_FromDouble(result);
  482. }
  483.  
  484. static int
  485. complex_nonzero(v)
  486.     PyComplexObject *v;
  487. {
  488.     return v->cval.real != 0.0 || v->cval.imag != 0.0;
  489. }
  490.  
  491. static int
  492. complex_coerce(pv, pw)
  493.     PyObject **pv;
  494.     PyObject **pw;
  495. {
  496.     Py_complex cval;
  497.     cval.imag = 0.;
  498.     if (PyInt_Check(*pw)) {
  499.         cval.real = (double)PyInt_AsLong(*pw);
  500.         *pw = PyComplex_FromCComplex(cval);
  501.         Py_INCREF(*pv);
  502.         return 0;
  503.     }
  504.     else if (PyLong_Check(*pw)) {
  505.         cval.real = PyLong_AsDouble(*pw);
  506.         *pw = PyComplex_FromCComplex(cval);
  507.         Py_INCREF(*pv);
  508.         return 0;
  509.     }
  510.     else if (PyFloat_Check(*pw)) {
  511.         cval.real = PyFloat_AsDouble(*pw);
  512.         *pw = PyComplex_FromCComplex(cval);
  513.         Py_INCREF(*pv);
  514.         return 0;
  515.     }
  516.     return 1; /* Can't do it */
  517. }
  518.  
  519. static PyObject *
  520. complex_int(v)
  521.     PyObject *v;
  522. {
  523.     PyErr_SetString(PyExc_TypeError,
  524.            "can't convert complex to int; use e.g. int(abs(z))");
  525.     return NULL;
  526. }
  527.  
  528. static PyObject *
  529. complex_long(v)
  530.     PyObject *v;
  531. {
  532.     PyErr_SetString(PyExc_TypeError,
  533.            "can't convert complex to long; use e.g. long(abs(z))");
  534.     return NULL;
  535. }
  536.  
  537. static PyObject *
  538. complex_float(v)
  539.     PyObject *v;
  540. {
  541.     PyErr_SetString(PyExc_TypeError,
  542.            "can't convert complex to float; use e.g. abs(z)");
  543.     return NULL;
  544. }
  545.  
  546. static PyObject *
  547. complex_conjugate(self, args)
  548.     PyObject *self;
  549.     PyObject *args;
  550. {
  551.     Py_complex c;
  552.     if (!PyArg_ParseTuple(args, ":conjugate"))
  553.         return NULL;
  554.     c = ((PyComplexObject *)self)->cval;
  555.     c.imag = -c.imag;
  556.     return PyComplex_FromCComplex(c);
  557. }
  558.  
  559. static PyMethodDef complex_methods[] = {
  560.     {"conjugate",    complex_conjugate,    1},
  561.     {NULL,        NULL}        /* sentinel */
  562. };
  563.  
  564.  
  565. static PyObject *
  566. complex_getattr(self, name)
  567.     PyComplexObject *self;
  568.     char *name;
  569. {
  570.     if (strcmp(name, "real") == 0)
  571.         return (PyObject *)PyFloat_FromDouble(self->cval.real);
  572.     else if (strcmp(name, "imag") == 0)
  573.         return (PyObject *)PyFloat_FromDouble(self->cval.imag);
  574.     else if (strcmp(name, "__members__") == 0)
  575.         return Py_BuildValue("[ss]", "imag", "real");
  576.     return Py_FindMethod(complex_methods, (PyObject *)self, name);
  577. }
  578.  
  579. static PyNumberMethods complex_as_number = {
  580.     (binaryfunc)complex_add, /*nb_add*/
  581.     (binaryfunc)complex_sub, /*nb_subtract*/
  582.     (binaryfunc)complex_mul, /*nb_multiply*/
  583.     (binaryfunc)complex_div, /*nb_divide*/
  584.     (binaryfunc)complex_remainder,    /*nb_remainder*/
  585.     (binaryfunc)complex_divmod,    /*nb_divmod*/
  586.     (ternaryfunc)complex_pow, /*nb_power*/
  587.     (unaryfunc)complex_neg, /*nb_negative*/
  588.     (unaryfunc)complex_pos, /*nb_positive*/
  589.     (unaryfunc)complex_abs, /*nb_absolute*/
  590.     (inquiry)complex_nonzero, /*nb_nonzero*/
  591.     0,        /*nb_invert*/
  592.     0,        /*nb_lshift*/
  593.     0,        /*nb_rshift*/
  594.     0,        /*nb_and*/
  595.     0,        /*nb_xor*/
  596.     0,        /*nb_or*/
  597.     (coercion)complex_coerce, /*nb_coerce*/
  598.     (unaryfunc)complex_int, /*nb_int*/
  599.     (unaryfunc)complex_long, /*nb_long*/
  600.     (unaryfunc)complex_float, /*nb_float*/
  601.     0,        /*nb_oct*/
  602.     0,        /*nb_hex*/
  603. };
  604.  
  605. PyTypeObject PyComplex_Type = {
  606.     PyObject_HEAD_INIT(&PyType_Type)
  607.     0,
  608.     "complex",
  609.     sizeof(PyComplexObject),
  610.     0,
  611.     (destructor)complex_dealloc,    /*tp_dealloc*/
  612.     (printfunc)complex_print,    /*tp_print*/
  613.     (getattrfunc)complex_getattr,    /*tp_getattr*/
  614.     0,                /*tp_setattr*/
  615.     (cmpfunc)complex_compare,    /*tp_compare*/
  616.     (reprfunc)complex_repr,        /*tp_repr*/
  617.     &complex_as_number,            /*tp_as_number*/
  618.     0,                /*tp_as_sequence*/
  619.     0,                /*tp_as_mapping*/
  620.     (hashfunc)complex_hash,     /*tp_hash*/
  621. };
  622.  
  623. #endif
  624.